#!/usr/bin/perl
 use strict;
# use warnings;
use Cwd;
use Data::Dumper;

#####################################################

#my ($USAGE) =  "\n\n\t****** USAGE of panam.pl PROGRAM ******\n\n\n\n\tUSAGE: perl $0 <panam.ini file> <taxonomic_group> \n\n\n\n";
#die "$USAGE" if ((scalar(@ARGV))< 1);
# my $option_file = $ARGV[0];
my $option_file = "panam.ini";
chomp($option_file);
my $taxo_group = $ARGV[0];
chomp($taxo_group);
die "\n\n\t=> Cannot find configuration file: $option_file.\n\n" unless (-e $option_file);
die "\n\n\t=> Configuration file, $option_file, appears to be empty!\n\n" if (-z $option_file);
open(OFILE, "<$option_file") || die "Cannot open input file : $option_file\n";

my $usearch;
my $NGS_id_Results;
my $panam_output = "panam_output";
# my $parsing;
my $query_seq; 
my $user_file;
my $preprocess_output = "preprocess_output";
my $seq_F;
my $seq_R;
my $primF;
my $primR;
my $trim =0;
# my $path = "bin/uclust3.0.617"; # Tung: no need
my $us_version = "3.0.617";
####
my $lib = "lib";
my $fast = "$lib/FastTree";
my $align = "$lib/hmmalign";
my $Reference = "Reference" ;
# my $build = "$lib/hmmbuild"; # Tung: no need
####
my $dom;
my %dom;
my $aff_norm;


while (<OFILE>) {
	chomp;
	my $line = $_;

# discard blank line
	if ($line =~ m/^\s*$/g) {
		next;
	}

# discard comment line
	elsif($line =~ m/^\s*#/g) {
		next;
	}
	
# get input options line
        else {
		
		if($line=~m/454_RUN_IDENTIFIER/gi) { 
			$line=~s/454_RUN_IDENTIFIER//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$NGS_id_Results = $line; 
			}
	 		die "\n\tOutput directory for NGS analyses is not defined. Check panam.ini.\n\n" unless ( $NGS_id_Results ne "");
			unless (-d $NGS_id_Results) { mkdir $NGS_id_Results || die "\n\tCould not create Directory $NGS_id_Results $!\n"; }
		}

# 		if (-d $NGS_id_Results."/".$panam_output) { `rm -r $NGS_id_Results/$panam_output` } 
		if (!(-d $NGS_id_Results."/".$panam_output)) { mkdir $NGS_id_Results."/".$panam_output || die "\n\tCould not create Directory $NGS_id_Results/$panam_output $!\n"; }


		if($line=~m/USER_FILE/gi) {
			$line=~s/USER_FILE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$user_file = $line; 
			}
			die "\n\n\t User file value missed. Check panam.ini.\n\n" unless (($user_file ne ""));
			die "\n\n\t Invalid user file value!\n\n" unless (($user_file eq "own") or ($user_file eq "pooledSamples") or ($user_file eq "eachSample"));	
		}


		if($line=~m/FORWARD_PRIMER_NAME/gi) {
			$line=~s/FORWARD_PRIMER_NAME//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$primF = $line; 
			}
		}
	
		if($line=~m/REVERSE_PRIMER_NAME/gi) {
			$line=~s/REVERSE_PRIMER_NAME//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$primR = $line; 
			}
		}

		if($line=~m/FORWARD_PRIMER_SEQUENCE/gi) {
			$line=~s/FORWARD_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$seq_F = $line; 
			}
		}
	
		if($line=~m/REVERSE_PRIMER_SEQUENCE/gi) {
			$line=~s/REVERSE_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$seq_R = $line;
			}
		}

		############### normalisation
		if($line=~m/NBR_SEQ_NORM/gi) {
			$line=~s/NBR_SEQ_NORM//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$aff_norm = $line; #print "aff --- $aff_norm\n";
			} 
		}
		
		#####################################

		if($line=~m/DOMAIN/gi) {
			$line=~s/DOMAIN//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) { 
				$dom = $line;			
				if (defined $dom ) { #push (@dom, $dom);
					$dom{$dom} = 1
				}
			}
			foreach my $d (keys %dom) {
				die "\n\n\t Domain name is not correct. Check panam.ini.\n\n" unless (($d eq "bacteria") or ($d eq "eukaryota") or ($d eq "archaea"));
			}	
		}
	}
}

die "\n\n\t Invalid number of sequences to pick for the normalization! Check panam.ini.\n\n" if ((defined $aff_norm ) and ($aff_norm !~ /^\d+$/));
if (!(defined $aff_norm)) {$aff_norm = '200' }


#print "$Reference\n";

# #############################################################
# 

############################################# trimmed profiles path ######################################################
# # # 
if (!(defined $primF) and (defined $primR))  {
	die "\n\n\t Primer forward not defined. Check panam.ini.\n\n";
}

if (defined $primF) {
	die "\n\n\t Primer reverse not defined. Check panam.ini.\n\n" unless (defined $primR);
	die "\n\n\t Primers' sequences missed. Check panam.ini.\n\n" unless ((defined $seq_F) and (defined $seq_R));
	die "\n\n\t Primers' sequences seem to be in the wrong format. Check panam.ini.\n\n" unless (($seq_F =~ /[ATCGUatcgu]/) and ($seq_R =~ /[ATCGUatcgu]/));
	$trim = 1
}

my @tax=();
my $Profiles;

if ($trim == 1) {
	#$Profiles = "Reference/Profiles_".$primF."_".$primR;
	$Profiles = "Reference/Trimmed_Profiles";
	if ((-d $Profiles)) {
		print "Trimmed profiles exist in $Profiles/\n";

	}
}

 
if ($trim == 0) {
	$Profiles= "Reference/Profiles";
}


# 
# ########################################################################################################
# 
# ########## hash reference sequences

unless (-e "$NGS_id_Results/$panam_output/Similarity_annotation/Sorted_Sequences") { `mkdir $NGS_id_Results/$panam_output/Similarity_annotation/Sorted_Sequences`; }


###################################

################################################################################### MODIF JCC 19/12/2011: 1 (sous)groupe taxonomique par noeud
#
#my $max=2;
#while (-e "$NGS_id_Results/$panam_output/Similarity_annotation/Sorted_Sequences/files_$taxo_group/file_cut_$max") {
#	$max++;	
#}
#
#$max--;
#my @num = (1 .. $max);
#
####################################################################################FIN MODIF JCC 19/12/2011 

# 
# # ##############################""
# 
if (!(-e "$NGS_id_Results/$panam_output/Alignment")) {
`mkdir $NGS_id_Results/$panam_output"/Alignment"`;
}


# MODIF JCC 19/12/2011 foreach my $fff (@num) {
#########Alignment

# MODIF JCC 19/12/2011 	if (-e "$NGS_id_Results/$panam_output/Similarity_annotation/Sorted_Sequences/files_$taxo_group/file_cut_$fff") { 
	if (-e "$NGS_id_Results/$panam_output/Similarity_annotation/Sorted_Sequences/files_$taxo_group/file_cut") { 
# MODIF JCC 19/12/2011 		my $ca = $align." -o ".$NGS_id_Results."/".$panam_output."/Alignment/".$taxo_group.".".$fff.".fasta --outformat a2m --mapali ".$Profiles."/fasta/".$taxo_group.".fasta ".$Profiles."/hmmprofil/".$taxo_group.".hmmprofil ". $NGS_id_Results."/".$panam_output."/Similarity_annotation/Sorted_Sequences/files_".$taxo_group."/file_cut_".$fff;	
		if (-e $NGS_id_Results."/".$panam_output."/Alignment/".$taxo_group.".fasta") {
			`rm $NGS_id_Results/$panam_output"/Alignment/"$taxo_group".fasta"`;
		}
		my $fastaProfiles = $Profiles."/fasta/".$taxo_group.".fasta";
		if (!(-e "$fastaProfiles")) {
			$fastaProfiles = $Profiles."/fasta/".substr($taxo_group, 0, -1).".fasta";
			if (!(-e "$fastaProfiles")) {
				$fastaProfiles = $Profiles."/fasta/".substr($taxo_group, 0, -2).".fasta";
			}
		}		
		my $hmmProfiles = $Profiles."/hmmprofil/".$taxo_group.".hmmprofil";
		if (!(-e "$hmmProfiles")) {
			$hmmProfiles = $Profiles."/hmmprofil/".substr($taxo_group, 0, -1).".hmmprofil";
			if (!(-e "$hmmProfiles")) {
				$hmmProfiles = $Profiles."/hmmprofil/".substr($taxo_group, 0, -2).".hmmprofil";
			}
		}		
		
		my $ca = $align." -o ".$NGS_id_Results."/".$panam_output."/Alignment/".$taxo_group.".fasta --outformat a2m --mapali $fastaProfiles $hmmProfiles ". $NGS_id_Results."/".$panam_output."/Similarity_annotation/Sorted_Sequences/files_".$taxo_group."/file_cut";	
		print "$ca\n";
		system "$ca 2>&1"  ;

# MODIF JCC 19/12/2011		print "Alignment files_$taxo_group.$fff done\n";
		print "Alignment files_$taxo_group done\n";
		
	}


# MODIF JCC 19/12/2011 }



